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1 Introduction 

In order to average or compare signals from 
functional brain images across subjects, e.g., in the 
case of fMRI, PET, or MEG, it is necessary to 
register the images with respect to each other. 
Registration techniques have been used also from 
TMS target area determination [1] to automatic 
segmentation of MRI images [2]. 

One of the advantages of spatially normalized 
images is that the locations of activation can be 
reported as Euclidean coordinates within a standard 
space. With segmentation, the transformation can be 
used to move the segments from one image to 
another: this way the geometrical knowledge of the 
anatomical information is automatically applied in 
the segmentation. 

Ultimately, inter-subject registration creates a 
correspondence map consisting of three parameters 
for each voxel in the image. These parameters 
represent the x, y, and z coordinates of the 
corresponding point in another image, typically a 
template image defined in the Talairach space. The 
number of free parameters would then be over 50 
million for an MRI image of size 256x256x256. 
Automated inter-subject registration procedures 
typically overcome the parameter estimation 
problem by parameterizing the deformation field in 
a lower dimension. Techniques that align the images 
using affine 1 transformations are widely used with 
brain atlases such as the one of Talaraich & 
Tournoux. The 12-parameter affine transformation 
can be determined either according to a set of 
control points given by the user, or automatically by 
using the AIR method as described in [3]. One 
method is to represent the transformation as a linear 
combination of smooth basis functions [4]. 

One important property of a useful spatial 
transformation is homeomorphicity, i.e., the 
transformation has to be one-to-one. The affine 
transformations are always homeomorphic, so that 
every point in the registered image has a unique 


1 An affine transformation is one that preserves straight 
and parallel lines. 


correspondence in the template image, and vice 
versa. The homeomorphism can not always be 
ensured with nonlinear transformations, except 
when the images are very close to each other [4]. 

We describe a method that can be used to determine 
a homeomorphic nonlinear transformation without 
floating point arithmetics. By using 8 and 16 bit 
integers throughout the program, the memory 
requirements can be decreased so that it is possible 
to determine a deformation with over 6 million 
parameters with a regular PC. The method is capable 
of handling large deformations between images, so 
it is not necessary to use an affine registration before 
calculating the nonlinear transformation matrices. 

2 Mathematical model 

The similarity measure S used is the difference of 
the gray-levels of the corresponding points in the 
images g x and g 2 , 

s = g 1 (x)-g 2 (x + t(x)) , (1) 

where t is the displacement vector defining the 
transformation. To simplify the equations and to 
ensure the symmetric processing of the images, the 
cost function E to be minimized is written as 


E = jjjhi (x + t(x))-g 2 (x-t(x))] 2 
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In addition to the similarity-dependent term, a 
smoothness constraint minimizing the derivatives of 
the displacement field is included in the cost 
function [5]. The smoothness required is controlled 
by parameter a. 


2.1 Minimization 

The necessary conditions for t to be an optimal 
solution can be found by the Euler-Lagrange 
equations [6], 
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where / is the integrand in Eq. (2), and the 
derivatives of the displacement field are denoted as 
(ti) x . = t( Idxj . From Eq. (2) we obtain 


= [gi(x+t)-g 2 (x-t)] • 
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Using Eqs. (4) and (5), Eq. (3) becomes 

[gl(x+t)-g 2 (x—t)] 

- a V 2 ?,(x) = 0 . 
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X + t(x) = x A + t A (x A ) 



Figure 1: Geometry of the transformation. 

The eigenvectors and eigenvalues satisfy also the 
equation 

m !=iw (10) 

and, on the other hand, according to the definition of 
matrix norm 

IIJtll=m f (11) 


Laplacian of the displacement field can be 
approximated with the term V 2 ^(x) ~ f*(x) - ^(x), 
where tfx) is the average of the values of the 
displacement component in the points around x 
[7]. The displacement field can then be calculated by 
an iterative process described by 


t l n+l = t l n + y i a(t~ti) 


[ gl (x+t n )-g 2 (x-t n )] 

dgls.,,^ , dg 2 

dxi 


,.(x + f) + ^(x-f) 


( 7 ) 


where t t n+l is the present, and t ” +1 is the new value of 
the component i of the displacement field; /is the 
relaxation parameter. 


2.2 Norm of the displacement field 

A mapping y = x ± t(x) is one-to-one, if and only if 
the determinant of the Jacobian of y is not zero, 
det(J y ) ^ 0. By denoting the Jacobian of the 
displacement field by J t , we obtain 


J y = I±Jt. (8) 

Let us show that the mapping y = x - t(x) is 
homeomorphic if the norm of the Jacobian of t is 
smaller than one, II Jt II < 1- The determinants of the 
Jacobians J t and J y are related by 

det(J y ) = - det(J t -1). (9) 

If det(J t - I) = 0, one of the eigenvalues of J t is 
equal to one (all satisfy the equation det(J t - A,*I) 
= 0 ). 


According to Eqs. (10) and (11), all eigenvalues are 
smaller than the norm. By limiting the norm to be 
smaller than one, we can ensure that det(J t - I) ^ 0, 
and the mapping is homeomorphic. 

The homeomorphism of the mapping y = x + t(x) 
can be shown in a similar manner. 

An additional advantage of limiting the norm of the 
Jacobian of the displacement field is, that the 
mapping x - t(x) —» x + t(x) can easily be converted 
into the form x A —» x A + t A (x) (see Fig. 1). Let us 
assume that we have a point x A on the image A (x A = 
x - t(x)), and we want to find the corresponding 
point on the image B (x B = x + t(x)). We can now 
search the point x by means of fixed-point iteration 

x" +1 = x A + t(x"), (12) 

where x" is the old and x" +1 is the new candidate for 
x. After solving x, we obtain the value of the 
mapping t A in the point x A directly from the equation 
t A (x A ) = 2t(x). 

2.3 Ensuring the homeomorphism 

It is clear that the smoothness constraint affects the 
norm of the Jacobian: when the derivatives of the 
displacement field are decreased, also II J t II 
decreases. Thus, we can ensure that the mapping 
stays homeomorphic by increasing the smoothness 
parameter a when II J t II grows close to one. 

The selection of the norm used does not affect the 
homeomorphism: 1-norm and infinity-norm are 
suitable choices because of the calculation 
simplicity. 




Figure 2: MRI images of two subjects. The crosses 
represent a pair of points that , according to the 
calculation result, correspond to each other. 


2.4 Implementation 

The algorithm was implemented as a collection of 
Matlab mex-files. Fixed-point arithmetics were 
used: the displacement field was determined as three 
matrices of signed 16-bit integers: 8 bits were used 
for the sign and the integer part, and the other 8 bits 
for the decimal part. 

Multi-resolution techniques were employed to speed 
up the processing. Each resolution step consists of 
the following operations: low-pass filtering of the 
MRI images; determination of the derivative 
matrices of the images (< dg/dx t in Eq. (7)); and 
iterating the displacements according to the Eq. (7). 
For MRI images of size 256x256x256, the basic 
form of the implementation requires memory for the 
three displacement fields of size 128x128x128 (16 
bits per element); two low-pass filtered images of 
size 256x256x256 (8 bits per element); and six 
derivative images of size 128x128x128 (8 bits per 
element). In addition, copies of the displacement 
fields are needed during the iterative process. Under 
Windows NT, 256 MB of memory ensures smooth 
operation even with on-line visualization of 
intermediate results. 

Current implementation uses linear interpolation in 
all resampling tasks. 

3 Evaluation 

The performance of the algorithm was evaluated by 
performing the registration between the MRI images 
of two subjects, and by creating and registering an 
artificial nonlinear transformation of an MRI image. 
The advantage of the latter evaluation method is that 
the correct answer is known precisely. 



Figure 3: Images before and after the registration 
and resampling. Registration error caused by the 
inhomogeneity can be seen in the throat area. 


3.1 Inter-subject case 

The algorithm was tested with a pair of MRI images 
of size 256x256x180. Mapping results for a dis¬ 
placement field of size 128x128x90 are visualized 
in Figs. 2 and 3. 

The calculation time in a inter-subject case is 10 to 
30 minutes, depending on the amount of dis¬ 
placement between corresponding objects. Pre¬ 
registration using a rigid body matching should be 
used to correct large rotational differences before the 
calculation of the displacement fields. This is partly 
because the orientation is not clearly visible in the 
low-pass-filtered images in the coarsest resolution, 
and partly because the norm of the Laplacian of the 
displacement field increases under rotations (II J t II = 
1 for a rotation of 90 degrees). 

The mapping result is most accurate in the areas 
with clear shapes and high gray-level gradients. In 
order to get a best possible mapping result, the 
histograms of the images should be adjusted so that 
the gray-level gradients on the edges of important 
objects would be clearly visible. 

3.2 Artificial case 

A nonlinear transformation of an MRI image was 
produced by creating a nonlinear displacement field 
using trigonometric functions (Fig. 4). Maximum 
displacement was 20 mm. 

Five control points were selected from the image, 
and the mapping error was tracked during the 
minimization procedure. The mapping error versus 
calculation time is visualized in Fig. 5. 






Figure 4: The artificial warping. The crosses mark 
a pair of points that correspond to each other 
according to the calculation result. 


4 Conclusion 

According to the performed tests, the algorithm 
seems to produce very accurate matching results 
even in the case of large deformations. 

In the current implementation, images are assumed 
to be of the same modality: measuring the similarity 
in terms of mutual information would make the 
algorithm more general [8]. 

The optimization was performed using Jacobi 
relaxation, i.e., all parameter values were updated 
simultaneously. Using Gauss-Seidel relaxation, i.e., 
updating the parameter values one by one, is in 
general more efficient than the Jacobi approach [9]: 
however, in this case the performance of Jacobi 
relaxation was clearly superior. Proper selection of 
the order in which the parameter values are updated 
could improve the results of the Gauss-Seidel 
relaxation. 
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Figure 5: Mapping error 
related to five control points as 
a function of the calculation 
time. Variance of the Gaussian 
filter and the distance between 
calculation grid points at each 
resolution level is marked on 
the top of the figure. Voxel size 
was Ixlxl mm 3 . 


































